7.3栅格数据可视化

7.3.1 环境矩阵数据可视化（二维）

##用x和y的边际概率密度函数绘制散点图：
# 同时生成三张图，包含散点图，对应x和y周到 边缘核密度图函数；，还可以加上对散点图的线性拟合；
## 参见网页：
https://github.com/ChrKoenig/R_marginal_plot

marginal_plot = function(x, y, group = NULL, data = NULL, lm_show = FALSE, lm_formula = y ~ x, bw = "nrd0", adjust = 1, alpha = 1, plot_legend = T, ...){
  require(scales)
  ###############
  # Plots a scatterplot with marginal probability density functions for x and y. 
  # Data may be grouped or ungrouped. 
  # For each group, a linear fit can be plotted. It is hidden by default, but can be shown by providing lm_show = TRUE.
  # The model can be modified using the 'lm_formula' argument.
  # The 'bw' and 'adjust' argument specify the granularity used for estimating probability density functions. See ?density for more information.
  # For large datasets, opacity may be decreased by setting alpha to a value between 0 and 1.
  # Additional graphical parameters are passed to the main plot, so you can customize axis labels, titles etc.
  ###############
  moreargs = eval(substitute(list(...)))

  # prepare consistent df
  if(missing(group)){
    if(missing(data)){
      if(length(x) != length(y)){stop("Length of arguments not equal")}
      data = data.frame(x = as.numeric(x), y = as.numeric(y))
    } else {
      data = data.frame(x = as.numeric(data[,deparse(substitute(x))]), 
                        y = as.numeric(data[,deparse(substitute(y))]))
    }
    if(sum(!complete.cases(data)) > 0){
      warning(sprintf("Removed %i rows with missing data", sum(!complete.cases(data))))
      data = data[complete.cases(data),]
    }
    group_colors = "black"
  } else {
    if(missing(data)){
      if(length(x) != length(y) | length(x) != length(group)){stop("Length of arguments not equal")}
      data = data.frame(x = as.numeric(x), y = as.numeric(y), group = as.factor(group))
    } else {
      data = data.frame(x = as.numeric(data[,deparse(substitute(x))]), 
                        y = as.numeric(data[,deparse(substitute(y))]),
                        group = as.factor(data[,deparse(substitute(group))]))
    }
    if(sum(!complete.cases(data)) > 0){
      warning(sprintf("Removed %i rows with missing data", sum(!complete.cases(data))))
      data = data[complete.cases(data),]
    }
    data = subset(data, group %in% names(which(table(data$group) > 5)))
    data$group = droplevels(data$group)
    group_colors = rainbow(length(unique(data$group)))
  } 

  # log-transform data (this is need for correct plotting of density functions)
  if(!is.null(moreargs$log)){
    if(!moreargs$log %in% c("y", "x", "yx", "xy")){
      warning("Ignoring invalid 'log' argument. Use 'y', 'x', 'yx' or 'xy.")
    } else {
      data = data[apply(data[unlist(strsplit(moreargs$log, ""))], 1, function(x) !any(x <= 0)), ]
      data[,unlist(strsplit(moreargs$log, ""))] = log10(data[,unlist(strsplit(moreargs$log, ""))])
    }
    moreargs$log = NULL # remove to prevent double logarithm when plotting
  }

  # Catch unwanted user inputs
  if(!is.null(moreargs$col)){moreargs$col = NULL}
  if(!is.null(moreargs$type)){moreargs$type = "p"}

  # get some default plotting arguments
  if(is.null(moreargs$xlim)){moreargs$xlim = range(data$x)} 
  if(is.null(moreargs$ylim)){moreargs$ylim = range(data$y)}
  if(is.null(moreargs$xlab)){moreargs$xlab = deparse(substitute(x))}
  if(is.null(moreargs$ylab)){moreargs$ylab = deparse(substitute(y))}
  if(is.null(moreargs$las)){moreargs$las = 1} 

  # plotting
  tryCatch(expr = {
    ifelse(!is.null(data$group), data_split <- split(data, data$group), data_split <- list(data))
    orig_par = par(no.readonly = T)
    par(mar = c(0.25,5,1,0))
    layout(matrix(1:4, nrow = 2, byrow = T), widths = c(10,3), heights = c(3,10))

    # upper density plot
    plot(NULL, type = "n", xlim = moreargs$xlim, ylab = "density",
         ylim = c(0, max(sapply(data_split, function(group_set) max(density(group_set$x, bw = bw)$y)))), main = NA, axes = F)
    axis(2, las = 1)
    mapply(function(group_set, group_color){lines(density(group_set$x, bw = bw, adjust = adjust), col = group_color, lwd = 2)}, data_split, group_colors)

    # legend
    par(mar = c(0.25,0.25,0,0))
    plot.new()
    if(!missing(group) & plot_legend){
      legend("center", levels(data$group), fill = group_colors, border = group_colors, bty = "n", title = deparse(substitute(group)), title.adj = 0.1)
    }

    # main plot
    par(mar = c(4,5,0,0))
    if(missing(group)){
      do.call(plot, c(list(x = quote(data$x), y = quote(data$y), col = quote(scales::alpha("black", alpha))), moreargs))
    } else {
      do.call(plot, c(list(x = quote(data$x), y = quote(data$y), col = quote(scales::alpha(group_colors[data$group], alpha))), moreargs))
    }
    axis(3, labels = F, tck = 0.01)
    axis(4, labels = F, tck = 0.01)
    box()

    if(lm_show == TRUE & !is.null(lm_formula)){
      mapply(function(group_set, group_color){
        lm_tmp = lm(lm_formula, data = group_set)
        x_coords = seq(min(group_set$x), max(group_set$x), length.out = 100)
        y_coords = predict(lm_tmp, newdata = data.frame(x = x_coords))
        lines(x = x_coords, y = y_coords, col = group_color, lwd = 2.5)
      }, data_split, rgb(t(ceiling(col2rgb(group_colors)*0.8)), maxColorValue = 255))
    }

    # right density plot
    par(mar = c(4,0.25,0,1))
    plot(NULL, type = "n", ylim = moreargs$ylim, xlim = c(0, max(sapply(data_split, function(group_set) max(density(group_set$y, bw = bw)$y)))), main = NA, axes = F, xlab = "density")
    mapply(function(group_set, group_color){lines(x = density(group_set$y, bw = bw, adjust = adjust)$y, y = density(group_set$y, bw = bw)$x, col = group_color, lwd = 2)}, data_split, group_colors)
    axis(1)
  }, finally = {
    par(orig_par)
  })
}

7.3.2 ggplot2可视化栅格

library(raster)
library(ggplot2)
library(rasterExtras)
library(RSpatial)
library(spocc)
library(dplyr)


occ <- scau

source("c:/Users/admin/Desktop/sdmenm.R")

bio1 <- raster("wc10/bio1.bil")

rang <- mcp(occ)



#plot
ggraster <- function(data,occ,mcp = FALSE){
  ## data为env形式：最好划定范围；
  ## occ仅包含经纬度信息；
  if(mcp == TRUE){
    rang <- mcp(occ)
    data <- raster::crop(data,rang)
  }else{
    data =data
  }
  data2 <- as.data.frame(data,xy=T)
  names(data2)[3] <- 'value'
  names(occ) <- c("longitude", "latitude")
  titll <- toupper(deparse(substitute(data)))

  ggplot() +
    geom_raster(data = data2, aes(x = x, y = y, fill = value)) +
    geom_point(data=occ, aes(x=longitude, y=latitude ),size= 5, col='green', cex=0.1) +
    coord_quickmap() +
    theme_bw() + 
    scale_fill_gradientn(colours=rev(terrain.colors(10)),
                         na.value = "#CCCCCC7f")+
    ggtitle(titll)
}
# ggraster(data,occ,mcp= TRUE)



## 第二种：对栅格数据分组后可视化：####
ggplotclass <- function(data){
  class <- summary(round(data,3)) %>% data.frame()
  ## 分别获取最小、第一四分卫，中值，第三四分卫，最大值；
  vacl <- class[,1][1:5] 
  ## 构造class:
  vacl2 <- c(vacl[1:2],1,vacl[2:3],2,vacl[3:4],3,vacl[4:5],4)

  ff <- reclassify(data2,vacl2)

  library(rasterVis)
  gplot(ff) +
    geom_raster(aes(fill = factor(value))) +
    coord_equal()+
    guides(color=guide_legend(title = deparse(substitute(data))))
}
# ggplotclass(data2)

7.3.3可视化物种概率密度池：

## 可视化物种概率密度池：
https://github.com/ChrKoenig/probpool/blob/master/demo/probpool.R

##
library(probpool)

### Create a new Probpool object from Ranunculaceae dataset
example_pp = probpool(env_pool = Ranunculaceae$Ran_env,
                          disp_pool = Ranunculaceae$Ran_disp,
                          occurrences = Ranunculaceae$Ran_occ,
                          interaction_matrix = Ranunculaceae$Ran_int, interaction_method = 1)

### Inspect Probpool
class(example_pp)
print(example_pp)
show(example_pp)
summary(example_pp)

### Plot Probpool
# 1. Plot entiree probpool
plot(example_pp, main = "Ranunculaceae dataset")

# 2. Plot individual species
plot(example_pp, focal_species = "Thal_simp", main = "Thalictrum simplex")

# 3. Plot individual spatial subset
plot(example_pp, focal_unit = c(1242,1241))

7.3.6 rastervis可视化栅格专包

## rastervis参见如下网站：
https://oscarperpinan.github.io/rastervis/

http://rstudio-pubs-static.s3.amazonaws.com/16756_56e9817f4cab42ceae5bc28b788c2c01.html

https://rpubs.com/alobo/rasterOnGM

## 基于高程线采用反距离加权法构建高程；
https://www.r-bloggers.com/2019/04/r-as-gis-for-ecologists/

https://www.rayshader.com/

7.3.7 可视化DEM

## 制作立体效果的2D地形图：
slope <- terrain(elevation, opt = "slope")
aspect <- terrain(elevation, opt = "aspect")
hill <- hillShade(slope, aspect, 40, 270)
plot(hill, col = grey(0:100/100), legend = FALSE, main = "Spain")
plot(elevation, col = rainbow(25, alpha = 0.35), add = TRUE)

7.7.8 可视化环境变量-maps

### BIO10 分类可视化 
bio10xx <- reclassify(bio10x,rcl=c(-Inf,24,1, 24,28,2, 28,30,3, 30,Inf,4))
m_class <- c(-Inf,24,1, 24,28,2, 28,30,3, 30,Inf,4)
bio10x %>% crop(.,range) %>% mask(.,range) %>% 
  reclassify(m_class) %>%  
  as.factor() -> r_class
levels(r_class)[[1]]$names <- c("24<","24-28","28-30",">30")
names(r_class) <- "BIO10"
r_class %>% 
  tm_shape() + 
  tm_raster(legend.show = TRUE,
            n = 4)+
  tm_legend(position=c("left","bottom"), frame=FALSE)+
  tm_shape(occ_data)+
  tm_symbols(col = "red", scale = .5,size=0.5,alpha = 0.7) +
  tm_shape(clip2)+
  tm_borders(col="gray50",alpha = 0.6)+
  tm_scale_bar(position=c("left", "top"),text.size = 0.8) + 
  #添加指北针
  tm_compass(type = "4star", position=c("right", "top")) +
  tm_xlab("Longitude",size = 1) +
  tm_ylab("Latitude",size = 1)+
  tm_graticules(ticks = TRUE,lines=FALSE)